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Abstract 

Google  search  algorithm  for  finding  relevant  information  on-line  based  on  keywords,  phrases, 
links,  and  webpages  is  analyzed  in  the  mathematical  and  sociological  settings  in  this  article.  We  shall 
first  survey  mathematical  study  related  to  the  Google  search  engine  and  then  present  a  new  analysis 
for  the  convergence  of  the  search  algorithm  and  a  new  update  scheme.  Next  based  on  sociological 
knowledge,  we  propose  to  use  in-  and  out-  linkages  as  well  as  use  the  second  order  linkages  to  refine 
and  improve  the  search  algorithm.  We  use  the  sociology  to  justify  our  proposed  improvements  and 
mathematically  prove  the  convergence  of  these  two  new  search  algorithms. 


1  Introduction 

Due  the  huge  volume  of  documents  available  through  the  world  wide  web,  a  search  engine  or  similar  ser¬ 
vice  is  now  absolutely  necessary  in  order  to  find  relevant  information  on-line.  About  12  years  ago,  Google, 
Yahoo  and  other  companies  started  providing  such  a  service.  It  is  interesting  to  know  the  mathematics  of 
computational  algorithms  behind  these  search  engines,  in  particular,  the  Google  search  engine  which  is 
so  successful  nowaday.  Many  mathematicians  and  computer  scientists,  G.  Golub  and  his  collaborators, 
Brezinski  and  his  collaborators,  Langville  and  Meyer,  among  other  pioneers  worked  on  the  mathematics 
of  the  Google  matrix  and  its  computation.  See  [25],  [16],  [9],  [8],  [31],  [32]  and  a  recent  survey  in  [2]. 
The  search  algorithm  also  becomes  a  very  useful  tool  in  many  Web  search  technologies  such  as  spam 
detection,  crawler  configuration,  trust  networks,  and  etc.,  and  find  many  applications  in  [12],  [13]  and 
etc..  Today,  it  is  the  most  used  computational  algorithm  in  the  world.  More  and  more  websites  such  as 
Facebook,  Amazon,  Netflix,  and  etc.,  use  similar  search  algorithms  for  various  purposes. 

As  almost  all  people  on  the  earth  use  the  Google  search  once  or  several  times  daily,  it  is  interesting 
to  see  if  the  search  results  make  any  sense  in  sociology.  How  could  it  be  adjusted  to  be  more  effective 
sociologically.  We  divide  this  article  into  two  parts.  We  shall  first  present  a  mathematical  description  of 
the  computational  algorithm  behind  the  Google  search  engine  and  summarize  the  recent  studies  related 
to  the  algorithm.  Then  we  present  a  sociological  analysis.  Based  on  the  sociology,  we  propose  two 
modifications  in  order  to  be  more  effective. 

This  research  is  supported  by  Army  Research  Office  Grant  #  W91  INF-1 1-1-0322. 

'mjlai@math.uga.edu.  Department  of  Mathematics,  The  University  of  Georgia,  Athens,  GA  30602. 

Lodawn@uga.edu.  Department  of  Sociology,  The  University  of  Georgia,  Athens,  GA  30602. 
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2  Mathematical  Analysis 


In  this  section  we  first  survey  the  mathematical  studies  related  to  the  Google  search  algorithm.  Then 
we  present  a  mathematical  analysis  which  is  new  to  the  best  of  our  knowledge.  Mainly  we  study  its 
convergence  as  it  is  an  iterative  method  since  the  Google  matrix  is  of  huge  size  so  that  any  direct  method  is 
impossible.  Then  we  shall  establish  the  convergence,  present  the  rate  of  convergence  and  a  local  updating 
scheme. 

2.1  Google  Search  Engine 

During  each  search,  after  entering  a  key  word  or  key  words,  the  Google  search  engine  will  produce  a 
ranking  vector  called  PageRank.  It  is  a  vector  of  huge  size  (several  billion  entries)  although  almost  all 
entries  arc  zero  in  general.  Each  entry  of  the  vector  presents  a  webpage.  Nonzero  value  in  an  entry  is  a 
weighting  factor  deciding  how  important  the  page  related  to  the  key  word(s)  just  entered.  The  importance 
is  determined  by  popularity.  The  ideas  for  this  computational  algorithm  have  been  reported  as  early  as 
1997-1998  by  Page  in  [37]  and  Brin  and  Page  in  [10].  See  also  [33]  for  explanation  of  the  ideas.  The 
algorithm  mainly  computes  a  numerical  weighting  (called  pagerank)  to  each  webpage  or  each  element 
in  a  hyperlinked  set  of  documents.  Pageranks  reflect  the  Google’s  view  of  the  importance  of  webages 
among  more  than  several  billion  webpages.  It  is  a  pragmatic  approach  to  improve  search  quality  by  going 
through  the  collective  intelligence  of  the  web  to  determine  a  page’s  importance. 

Let  v  be  a  vector  of  with  N  >  8  billion.  Any  unit  vector  in  MA  is  associated  with  a  webpage. 
Mathematically,  Pagerank  is  a  probability  function  of  vectors.  For  a  page  v,  let  P(v)  be  the  PageRank  of 
v.  By  a  stroke  of  the  ingenious,  Brin  and  Page  thought  that  P(v)  is  the  weighted  sum  of  the  PageRanks 
of  all  pages  pointing  to  v.  That  is, 

P(v)=  ^P(u)((u,v),  (1) 

u£ytv 

where  Av  is  the  collection  of  all  pages  connected  to  v  and  /:'(u.  v)  e  [0, 1]  be  a  weight  dependent  on  the 
number  of  links  from  page  u  to  the  other  pages.  More  precisely,  let  L(u)  be  the  number  of  linkages  from 
u  and  then  define  t(u.  v)  =  1  / L(u)  if  u  6  Av  and  0  otherwise  for  any  other  page  u.  Such  a  definition 
is  cyclic.  Also,  a  web  surfer  may  get  bored  and  abandons  the  current  search  by  teleporting  to  a  new  page. 
By  a  second  stroke  of  ingenious,  Brin  and  Page  introduced  a  damping  factor  and  modified  the  Pagerank 
function  P(v)  to  be 

P(v)  =  d  ^2  P(u)l(u,  v)  +  1  Vunit  vector  v  €  RN ,  (2) 

ueAv 

where  d  €  (0, 1)  is  a  damping  factor.  It  indeed  leads  to  a  convergent  algorithm  and  the  PageRank  P  will 
be  unique  defined.  In  summary,  we  have 

•  PageRank 

P  =  [P(Vl),--- ,P(V7V)]T.  (3) 

•  M  =  [((v,.  Vj  )]  1<);  .<N  be  the  adjacency  matrix  of  size  N  x  N  with  convention  i(vt.  v4)  =  0. 

•  1  =  [1, 1,  -  -  -  ,1]T 
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•  The  Google  Matrix 


5  =  dM  +  ^j^  1  1T.  (4) 

which  is  a  stochastic  matrix,  where  d  >  0  is  a  damping  factor.  In  Google  search  engine,  d  =  0.85. 

To  compute  the  PageRank  P,  the  Google  search  engine  computes:  starting  with  an  inital  vector 
Pi,e.p.,Pi  =  1/N, 

Pjfe+1  =  0Pfc,Vfc  =  l,2,---  (5) 

until  the  difference  of  P^+i  and  P&  within  a  given  tolerance. 

For  an  elementary  explanation  of  Google  computational  algorithm,  see  a  monograph  [33].  Many 
mathematical  problems  arise  from  the  algorithm.  For  example,  does  the  computation  (5)  converge  and 
what  is  the  convergence  rate?  Is  the  computation  stable?  How  to  speed  up  the  iterations  in  (5)?  How  to 
do  the  computation  in  parallel?  Can  we  have  d  =  1  ?  How  can  we  add  or  delete  a  page  web  from  the 
PageRank  function? 

2.2  Summary  of  Mathematical  Results 

The  major  question  for  the  computational  algorithm  (5)  in  the  previous  subsection  is  the  convergence  rate 
of  the  iterations  in  (5).  The  iteration  is  a  power  method.  The  convergence  is  dependent  on  the  second 
largest  eigenvalue  as  the  first  largest  eigenvalue  is  1  (cf.  [27]).  The  value  of  the  second  largest  eigenvalue 
or  eigenvalues  is  the  damping  factor  as  shown  first  in  [24] .  This  fact  is  later  established  in  [39]  by  using 
Jordan  canonical  form  of  the  Google  matrix.  It  was  also  proved  in  [14]  where  the  researcher  showed 
that  the  second  largest  eigenvalue  is  the  dampling  factor  even  the  Google  matrix  has  more  than  1  largest 
eigenvalues.  It  is  also  proved  in  [30].  In  [1 1],  the  searcher  interpreted  the  PageRank  algorithm  as  a  power 
series  which  leads  to  a  new  proof  of  the  several  known  results  about  the  PageRank,  e.g..,  the  convergence 
rate  is  the  exactly  the  damping  factor  by  [24]  and  presented  a  new  inside  about  PageRank  which  is  the 
probability  that  a  specific  type  of  random  surfer  ends  his  walk  at  the  page.  Usually,  the  random  surfer  is 
assumed  to  be  uniform  distribution  suggested  by  Page  and  Brin.  In  fact,  all  the  properties  of  the  PageRank 
hold  for  arbitrary  probability  setting.  We  shall  present  another  proof  in  the  next  subsection 

A  second  question  is  if  the  computation  is  stable  or  what  the  condition  number  of  the  Google  matrix 
is.  In  [18],  the  researchers  viewed  the  PageRank  vector  to  the  stationary  distribution  of  a  stochastic  matrix, 
the  Google  matrix.  They  analyzed  the  sensitivity  of  PageRank  to  changes  in  the  Google  matrix,  including 
addition  and  deletion  of  links  in  the  web  graph  and  presented  error  bounds  for  the  iterates  of  the  power 
method  and  for  their  residuals.  See  also  various  properties  in  [9]  including  the  condition  number  of  the 
Google  matrix. 

A  central  question  is  how  to  accelerate  the  Google  computational  algorithm.  There  are  several  ap¬ 
proaches. 

•  Extrapolation  In  [8],  the  researchers  studied  the  computation  of  the  nonnegative  left  eigenvector 
associated  with  the  dominant  eigevalue  of  the  PageRank  matrix  and  proposed  to  use  extrapolation 
methods  to  approximate  the  eigenvector.  Furthermore,  in  [9],  the  researchers  proposed  some  ac¬ 
celeration  methods  due  to  the  slow  convergence  of  the  Power  method  for  the  PageRank  vector. 
The  acceleration  metod  can  be  explained  as  the  method  of  Vorobyev  moments  and  Krylov  sub¬ 
space  method.  In  addition,  many  properties  about  the  Google  matrix  and  PageRank  vectors  were 
presented. 
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•  Sparsity  Split  In  [40],  these  Chinese  researchers  use  the  sparsity  of  the  hyperlink  matrix  (one  of  the 
two  matrices  in  Google  Matrix)  and  proposed  a  linear  extrapolation  to  optimize  each  iteration  in 
the  Power  method  for  the  computation  of  PageRank  and  reduction  of  the  storage  space.  The  idea  is 
simple  and  the  engineers  in  Google  may  have  already  implemented. 

•  Iterative  Solutions  The  iteration  in  (5)  can  be  reinterpreted  as  a  linear  system  to  be  given  in  (8).  As 
explained  in  [15],  the  Google  matrix  is  very  large,  sparse  and  non- symmetric.  Solution  of  the  linear 
system  by  a  direct  method  is  not  feasible  due  to  the  matrix  size.  Many  iterative  solutions  such  as 
Gauss-Jacobi  iterations.  Generalize  Minimum  Residual  (GMRES),  Biconjugate  Gradient  (BiCG), 
Quasi-Minimal  Residual  (QMR),  Conjugate  Gradient  Squared  (CGS),  Biconjugate  Gradient  Sta¬ 
bilized  (BiCGSTAB),  Chebyshev  Iterations.  The  researchers  proposed  to  use  the  parallel  Block 
Jabobi  iteration  with  adaptive  Schwarz  preconditioners. 

•  Parallel  Iteration  In  [15],  the  researchers  described  a  parallel  implementation  for  PageRank  com¬ 
putation  and  demonstrated  that  the  parallel  PageRank  computing  is  faster  and  less  sensitive  to  the 
changes  in  teleportation. 

•  Lumping  Dangling  Nodes  In  [19],  the  researchers  lumped  all  of  the  dangling  nodes  together  into 
a  single  node  and  the  Google  matrix  can  be  reduced.  They  showed  that  the  reduced  stochastic 
matrix  has  the  same  nonzero  eigenvalues  as  the  full  Google  matrix  and  the  convergence  rate  is  the 
same  when  the  Power  method  is  applied.  In  [34],  the  researchers  further  reduced  the  Google  matrix 
by  lumping  more  nodes  which  arc  so-called  weakly  nondangling  nodes  together  and  the  reduced 
matrix  has  the  same  nonzero  eigenvalues  as  the  Google  matrix. 

•  Distributed  Randomized  Algorithm  In  a  series  of  papers  ([20],  [21],  [22],  the  researchers  applied 
distributed  randomized  algorithm  for  PageRank.  The  idea  is  to  let  each  page  compute  its  own 
value  by  exchanging  information  with  its  linked  pages.  The  communication  among  the  pages  is 
randomized  so  as  to  make  it  asynchronous.  The  communication  among  pages  in  the  sense  that  data 
transmitted  over  the  links  is  received  correctly  or  could  have  some  noises. 

There  are  several  other  improvements.  For  example. 

One  is  to  update  the  PageRank.  In  [33],  the  researchers  considered  the  PageRank  vector  to  be  an 
state  of  homogeneous  irreducible  Markov  chain  and  the  chain  requires  updating  by  altering  some  of  its 
transition  probabilities  or  by  adding  or  deleting  some  states.  They  proposed  a  general  purpose  algorithm 
which  simultaneously  deals  with  both  kinds  of  updating  problems  and  proved  its  convergence.  We  shall 
present  another  approach  for  local  updating  to  be  discussed  in  the  next  subsection. 

Another  improvement  is  to  distinct  the  webpage  with  equal  rankings.  In  [35],  the  researchers  proposed 
an  idea  to  organize  the  search  result  produced  by  the  PageRank  where  several  pages  have  the  same  rank 
score  so  that  the  surfer  can  get  more  relevant  and  important  results  easily.  Their  idea  is  to  add  a  weight  to 
each  page. 

We  finally  remark  that  besides  the  PageRank  algorithm,  there  arc  many  other  search  algorithms  avail¬ 
able.  For  example,  Kleinberg’s  HITS  algorithm  is  a  classic  Hypertext  Induced  Topics  Search  (HITS) 
algorithm  invented  by  Kleinberg  (cf.  [28]).  Kleinberg  viewed  that  most  pages  arc  a  hub  and  authority, 
i.e.,  a  general  page  not  only  points  to  other  pages  and  also  receives  points  from  other  pages.  Again,  the 
HITS  algorithm  uses  the  power  method  to  the  hub  weight  vector  and  authority  weigh  vector.  See  [33]  and 
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[6]  for  more  explanation.  For  another  example,  Yahoo  uses  another  type  of  algorithms.  See  last  section 
of  this  paper. 

In  summary,  the  study  of  Google  computational  algorithm  requires  many  mathematical  tools  such 
that  graph  theory  (cf.  [?]),  random  walks  and  Markov  chains  (c.g.,  [4],  [33]),  numerical  analysis  (cf. 
[16]),  [15],  [8]),  stochastic  matrices  (cf.  [18]),  networks  theory  (cf.  [6]),  etc..  See  [2]  for  a  survey 
of  recent  techniques  for  on  link-based  ranking  methods  and  the  PageRank  algorithm  to  find  the  most 
relevant  documents  from  the  WWW. 

2.3  Our  Convergence  Analysis 

We  see  that  P(v)  >  0.  We  now  claim  that  P(v)  <  1.  In  fact,  we  can  prove  the  following 

Theorem  2.1  Let  v*,  i  =  1,  •  •  •  ,N  be  the  standard  unit  vectors  in  M,v.  Then  the  PageRanks  P(v,;) ,  i  = 
1,  •  •  •  ,  N  satisfy 

N 

Yp^)  =  1- 

i= 1 

Proof.  Indeed,  as  in  the  definition  (2),  P(vj),  i  =  1,  •  •  •  ,  N  arc  dependent  on  each  other.  We  recall 

P  =  [-P(vr),  -  -  -  ,P(wv)]T 


is  the  PageRank  vector  in  M;V .  We  can  see  that 

N 

2^(vj,v*)  =  1  (6) 

i= 1 

for  j  =  1,  •  •  •  ,  N.  We  call  £(vi,  v ; )  adjacency  functions.  It  follows  from  (2),  we  can  write 

P  =  —jy—  1  +  dMP ,  (7) 

where  M  =  [((v,,  v?)]  |<?;  -<N  be  the  adjacency  matrix  of  size  N  x  N  with  convention  ('(v,.  Vj)  =  0  and 
1  =  [1, 1,  •  •  •  ,  1]T  €  Hence,  we  have 

N  , 

P(vi)  =  1TP  =  +  dlTMP 

1=1 

N 

=  1  —  d  +  dlTP  =  1  —  d  +  d'Y,  P(vi), 

1=1 


N 

where  we  have  used  (6)  to  have  1 T M  =  1T.  That  is,  we  have  (1  —  d)  Yj  P(v*)  =  1  —  d  and  hence, 

i= 1 
N 

Y^  P(vi)  =  1  since  d  f  1.  This  completes  the  proof.  ■ 

i=  1 
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The  PageRank  algorithm  offers  two  ways  to  calculate  the  PageRank  vector  P:  one  is  by  algebraic 
method  and  the  other  by  an  iterative  method.  Indeed,  by  (7),  we  have 

(I-dM)  P  =  ^l  or  P=1^d(I-dM)~1l,  (8) 

where  I  is  the  identity  matrix  of  size  N  x  N.  This  algebraic  method  will  be  well-defined  if  I  —  dM  is 
invertible.  Indeed,  we  have 

Theorem  2.2  Let  I  be  the  identity  matrix  and  M  be  the  adjacency  matrix  defined  above.  If  d  €  (0, 1), 
then 

I-dM 


is  invertible. 


Proof.  We  first  recall  a  matrix  A  =  «y]  i<ij<Ar 's  diagonally  dominant  if 


N 


3  =  1 


for  alH  =  1,  •  •  •  ,  N.  It  is  easy  to  see  that  (I  —  dM)  is  diagonally  dominant  matrix  by  using  (6)  with 
d  <  1.  It  is  known  (cf.  [27],  p.  178)  that  every  diagonally  dominant  matrix  is  nonsingular,  and  hence  is 
invertible.  So  is  I  —  dM.  m 

Certainly,  finding  the  inverse  ( I  —  dM)~l  is  expensive  in  computation.  An  easy  approach  is  to  use  an 
iterative  method,  for  example,  power  method.  Starting  with  Po  =  1  /N,  we  iteratively  compute 

Pfc  =  GPk-1  =  '  v  f/ •  +  dMPfc_!  (9) 

by  using  Theorem  2.1  for  k  =  1, 2,  •  ■  •  We  now  show  that  P / ,  k  >  1  converge.  Let  us  use  the  l\  norm 
for  matrices  to  measure  the  errors  P&  —  P.  Recall  the  t\  norm  for  matrix  A  with  A  =  [ciij]i<i,j<N  is 


N 

i  =  max  >  |  a. 

7=i>"'  i  tv 

i=i 


%j\- 


Then  we  have 


Theorem  2.3  Let  P  be  the  PageRank  vector  defined  by  the  algebraic  method.  For  all  k  >  1, 


|Pfe  -  P||i  <Cd\ 


where  C  is  a  positive  constant. 

Proof.  It  is  easy  to  see  from  (7)  and  the  definition  (9)  of  iterations  that 

Pfc  -  p  =  dM( ;pfc_!  -  P)  =  •  •  •  =  (dM)k( P0  -  P). 
Thus,  letting  C  =  ||Po  —  P||i, 

||Pfc  -  P||i  <  dk\\Mk\hC  <  Cdk\\M\\i. 
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Finally  we  note  that  ||M||i  =  1  by  using  (6)  to  complete  the  proof.  ■ 

We  remark  that  the  convergence  of  P/  to  P  was  shown  numerically  in  [37]  for  full  size  and  a  half 
size  link  databases.  The  rates  of  convergence  arc  nearly  the  same  for  these  two  sizes.  Also  in  [23], 
the  contributor(s)  of  this  wiki  page  uses  the  power  method  for  the  leading  eigenvector  to  explain  the 
convergence  of  the  iterative  method  above  without  giving  a  convergence  rate  as  it  is  based  on  the  ratio 
of  the  second  largest  and  the  largest  eigenvalues  of  the  Google  matrix.  The  second  eigenvalue  of  the 
Google  matrix  was  studied  in  [24]  where  the  rate  of  convergence  of  PageRank  was  determined  to  be  0.85, 
the  damping  factor  d.  As  mentioned  in  a  previous  section,  the  convergence  rate  has  been  determined  by 
several  methods  already.  Our  result  in  Theorem  2.3  gave  another  proof.  It  is  much  simpler  and  easier 
fashion  than  the  study  in  [24]. 


2.4  A  New  Local  Update  Scheme 

We  discuss  how  to  do  local  updates.  That  is,  in  practice,  webpages  arc  updated  constantly  in  the  sense 
that  the  number  L(vj)  of  links  from  webpage  vt  changes  very  often,  usually  increases.  A  new  adjacency 
matrix  M  is  needed  updated  everyday.  Instead  of  using  an  updated  matrix  M  to  compute  a  new  PageRank 
vector  P  for  all  components,  we  discuss  an  approach  to  update  partial  components  of  P.  We  need  the 
following  well-known  formula: 


Lemma  2.1  (Sherman-Morrison’s  formula, 1949)  If  A  is  invertible  and  x,  y  are  two  vectors  such  that 
yTA~1x  7^—1,  then  A  +  xyT  is  invertible  and 


(A  +  xyT)~1  =  A-1 


A-1xyTA_1 
1  +  yTA-1x' 


Assume  that  we  have  an  updated  £4  (vz,  v7)  for  a  fixed  integer  j  and  M+  be  the  updated  adjacency 
matrix.  Let  n  =  n{j)  and  n+  =  n+  (j )  be  the  previous  and  updated  numbers  L(\j )  of  links  from  webpage 
Vj.  Let  us  forcus  the  case  that  only  more  webpage  links  are  added  to  the  existing  links  from  the  webpage 
Vj.  In  this  case  n+  >  n.  Hence,  we  have 

M+  =  M  +  XjvJ, 


where  Xj  =  [£+(vi,Vj)  —  £(vi,vj),i  =  1,  •  •  •  ,  N]T  with  £+(vi,  Vj)  is  the  updated  version  of  £(vi,Vj) 
and 

,, , ,  „ ,  ..  ,  n_i_  —  n  n+  —  n 

M  —  M+  1  =  1  —  n  n. i_  4 - <  2 - . 

n+  n+ 

We  shall  use  P  +  to  be  the  updated  PageRank  vector  satisfying 


P+ 


1  -d 
N 


1  +  dM+  P 


By  Theorem  2.2,  we  know  I  —  dM+  is  invertible.  By  the  inverse  formula  in  Lemma  2.1,  we  know 
vj  (I  —  dM)~1xj  f  —1  and  hence,  letting  a  =  1  +  v  J  (I  —  dM)~1xj  f  0, 

(1  -  dAf+r1  =  (1  -  dM)~l  -{I-  dM)~1xjvJ  (1  -  dM)~l/a.  (10) 
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It  follows 


P+  =  (l-dM+)~1^l 

=  (1  -dMy1^-^!-  i(/-dM)_1Xjv7(l-dM)-1^4l 

=  P-  -(/-dM)-1x7-vT P 
a  J 

=  P  -  -^(J-  dM)-1*-.  (11) 

a 

p. 

where  Pj  is  the  jth  component  of  vector  P.  Let  us  write  e  =  —  (I  —  dM  'f  'x;  for  the  update  for  P  to 

obtain  P  and  we  next  discuss  how  to  find  the  update  efficiently  as  solving  (/  —  dM)~  'xy  is  expensive 
in  computation. 

We  need  the  following  (cf.  [27],  p.  198) 

Lemma  2.2  (Neumann  Series)  If  A  is  an  n  x  n  matrix  and  ||A||i  <  1,  then  I  —  A  is  invertible  and 

OO 

a-Ar^y^AK 

k= 0 

As  shown  before  ||M||i  =  1  and  d  G  (0, 1),  letting  A  =  dM  in  Lemma  2.2  above,  we  have 

OO 

(/  -  dM)-lKj  =  dkM%. 

k= 0 

Algorithm  2.1  (A  Single  Update)  Choose  an  integer  K  large  enough  such  that  dK  is  within  a  given 
tolerance.  Start  with  y  =  x;  and  for  k  =  1,  •  •  •  ,  K  —  1,  we  compute 

y  =  dMy  +  xj. 

Then  writing  y  =  [yi,  ■  ■  ■  ,  y_/v]T»  we  let  a  =  1  +  yj  and  e  =  (Pj / ct)y.  Thus,  by  (11),  P+  «  P  —  e  within 
the  given  tolerance. 

The  above  algorithm  leads  to  sequentially  updating  P,  i.e.  update  one  webpage  or  one  keyword  or 
one  document  at  a  time.  See  Remark  5.1  for  a  parallel  update  scheme. 

3  Sociological  Analysis 

Sociologists  have  studied  the  linkages,  networks  for  a  long  time.  Many  concepts  of  centrality  measures 
have  been  created.  In  [7],  Bonscich  proposed  a  measure  of  centrality  c(a,  3)  with  normalization  factor  a 
and  parameter  fl  which  is  a  function  of  unit  vectors  of  R;V.  For  each  unit  vector,  say,  v*, 

Ci(a,3 )  =  +  3cj(a,3))Rij 

j 
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(12) 


where  Rij  is  an  entry  of  matrix  relation  R  which  may  not  be  symmetric  and  main  diagonal  elements  arc 
zero.  In  matrix  notation 

c(a,/3)  =  a(I  —  /3R)_1R1.  (13) 

See  page  1173  of  [7].  If  we  use  the  adjacent  matrix  for  R  and  normalization  factor  a  =  (1  —  d)/N ,  we 
see  that  (12)  is  just  (2)  when  ft  =  d,  the  damping  factor.  If  a  =  (1  —  d)/N  and  R1  =  1,  we  can  see  that 
(13)  is  the  same  as  the  equation  on  the  right-hand  side  of  (8).  Using  these  two  parameters,  the  PageRank 
is  just  c(a,  ft).  According  to  Bonacich,  when  R  is  asymmetric,  e.g.,  R  =  M ,  c (a,  (5)  measures  prestige 
if  ft  >  0.  This  gives  a  sociological  justification  that  the  Google’s  PageRank  does  measure  the  important 
relevance  scores  for  each  webpage. 

Note  that  when  R  is  symmetric,  c(a..  ft)  measures  centrality.  Also,  as  in  [7],  ft  can  be  negative.  For 
example,  in  a  communication  network. 

As  pointed  out  by  Bonacich,  his  measure  c (cr,  ft)  seems  hopelessly  ambiguous,, 


4  Improved  Search  Algorithms 

We  now  explain  the  improvements  discussed  in  the  previous  section  in  mathematical  terminology.  Recall 
all  the  notations  and  definitions  in  a  previous  section.  We  use  page  to  denote  a  webpage,  a  keyword  or 
phrase,  a  document,  etc.. 


4.1  Convergence  Analysis  of  In-  and  Out-Linkage  Search  Algorithm 

For  each  page  v,  Av  is  the  collection  of  all  pages  connected  from  v.  We  now  search  all  u  £  Av  to  see  if 
u  connects  to  v  also.  We  denote  by  Aft  the  collection  of  such  page  u  £  Av.  Note  that  Aft  could  be  an 
empty  set  if  v  is  an  unknown  page  to  the  world.  Anyway,  Aft  C  Av  is  a  subset.  Let 

L(v)  =  #(Av)  +  #(A+),  (14) 


where  #(A)  stands  for  the  number  of  entries  in  set  A.  It  is  called  the  cardinality  of  A.  We  define  the  new 
adjacency  entry 


*(v,u) 


'VUyft 
<  2/L(v), 
0 


if  u  is  not  connected  to  v  from  u 
if  u  is  connected  to  v 
otherwise. 


(15) 


According  to  our  Plan  1,  the  new  PageRank  function  P(v)  is  now  defined  by 


P(v)  =  i  Y.  f(u«".')+X' 

usAv 


V  any  unit  vector  v  £  M;Y, 


(16) 


where  d  £  (0, 1)  is  a  damping  factor.  To  find  the  value  of  P(v),  we  use  similar  iterative  algorithm  as  in 
Section  2.  Let 


M  = 


-I 


(17) 


We  are  now  ready  to  define 
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Algorithm  4.1  (In-  and  Out-linkage  Search  Algorithm)  Starting  from  Po  =  1  /N,  we  iteratively  com¬ 
pute 

Pk  =  ^-^1  +  dMPk-u  Vfc  =  1,2,, ....  (18) 

For  convenience,  we  write  P  to  the  new  PageRank  vector  of  all  pages  Vj,  <E  R;V,  i  =  1,  •  •  •  .  N . 

Based  on  the  analysis  in  Section  2,  we  can  conclude  the  following 

Theorem  4.1  Let  P  be  the  PageRank  vector  defined  above.  For  all  k  >  1, 

||P*  -P||i  <  Cdk,  Vk>  1, 

where  C  is  a  positive  constant  independent  ofk. 


4.2  Convergence  Analysis  of  Second  Order  linkage  Search  Algorithm 

Next  we  study  our  improved  search  algorithm  2.  In  addition  to  use  Af,  we  count  the  links  from  u  £  Av. 
That  is,  we  look  at  Au.  More  number  of  links  in  ,4U.  the  less  contribution  to  u  according  to  the  sociological 
justitication  in  the  previous  section.  Thus,  we  define 


n(v) 


#(av)  +  y 

usAv 


1 

#(A0 


and 


?(u,v) 


1  _i _ 1_ 

1  +  #(AU) 
*  n(v) 

0 


if  u  £  Av 
otherwise. 


It  is  easy  to  see 

Lemma  4.1  For  any  page  v 


N 

£?(v;>v)  =  l. 


3=1 

Thus  the  new  PageRank  function  P(v)  is  now  defined  by 


(19) 


(20) 


(21) 


p(v)  =  d  y  p(u)?(U,v)  +  i^, 

usAv 


V  any  unit  vector  v  £  ~RN , 


(22) 


where  d  £  (0, 1)  is  a  damping  factor.  To  find  the  value  of  P(v),  we  do  the  iterations  as  before.  Letting 


M  = 


Z(vi,Vj 


(23) 


Algorithm  4.2  (Second  Order  linkage  Search  Algorithm)  Starting  from  Po 
compute 


Pk 


1  -d 
N 


1  +  dMPfc-i , 


Vk  =  1,2,,.... 


1/N,  we  iteratively 
(24) 
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Similarly,  we  write  P  to  the  new  PageRank  vector  of  all  pages  v(,  e  M.N ,  i  =  1,  ■  •  •  ,  N. 

Based  on  the  analysis  in  Section  2,  we  can  conclude  the  following 

Theorem  4.2  Let  P  be  the  PageRank  vector  defined  above.  For  all  k  >  1 , 

||Pfc  -P||i  <  Cdk,  \/k>  1, 

where  C  is  a  positive  constant  independent  of  k. 

5  Yahoo  Search  Engine 

Yahoo  search  engine  is  quite  different  from  the  Google  search  engine  conceptually  and  hence,  computa¬ 
tionally.  The  following  description  of  the  Yahoo  search  engine  is  based  on  [38]. 

Yahoo  emphasizes  on  the  similarity  of  any  two  webpages,  two  documents,  or  two  web  hubs.  There  are 
many  sources  of  similarity  information  including  link  similarity,  text  similarity,  multimedia  component 
similarity,  click  through  similarity,  title  similarity,  URL  similarity,  cache  log  similarity,  and  etc..  For 
simplicity,  let  us  consider  the  search  engine  for  documents.  Even  for  documents,  there  are  many  categories 
of  similarity  such  as  title  similarity,  author  similarity,  text  similarity  and  etc..  For  example,  the  text 
similarity  is  the  dot  product  of  two  vectors  which  are  the  frequency  usages  of  various  words.  That  is, 
suppose  that  two  documents  use  total  M  words.  For  one  document,  let  x  be  the  vector  of  length  M  whose 
each  entry  is  the  frequency  usage  of  a  word  in  the  document.  Similar  for  vector  y  for  the  other  document. 
Then  x  •  y  denote  the  similarity  between  these  two  documents. 

Suppose  that  the  base  set  of  all  on-line  documents  contains  N  documents.  The  global  similarity 
matrix  W  =  ['Wifi  i<i,j<N  with  Wij  >  0  representing  the  strength  of  the  similarity  between  document  i 
and  document  j.  Assume  that  v'tJ  <  M  for  all  i,j  =  1,2,---  ,  N.  Suppose  that  there  arc  K  categories 
of  various  similarities.  For  each  k  =  1,  •  •  ■  ,  K,  let  x(fc)  =  [aq*.,  aq k,  •  ■  ■  ,  x^k\  F  be  a  vector  with 
nonnegative  entries  satisfying 

N 

^2 Xik  =  l. 

i= 1 

It  is  a  vector  of  scores  or  confidence  values  for  all  documents  with  respect  to  the  klh  category  of  similarity. 
Note  that  the  global  similarity  matrix  combines  all  sources  of  similarity  information  together.  Thus, 
=  f(wij(  1),  •  •  •  ,  Wij(K ))  with  Wij(k),  k  =  1,  •  ■  ■  ,  K  being  the  similarity  between  document  i  and 
document  j  with  respect  to  the  category  k.  Yahoo  considers  a  global  similarity  objective  function 

K  N 

P(x)  =  £  x(fc)TIUx(fc)  =  EE  Xjk  liUj  X jk  • 

k= 1  k= 1  i,]= 1 

Yahoo  engine  will  have  to  first  do  data  preparation  or  preprocessing  by  generating  the  global  similarity 
matrices  W  based  on  training  and  learning  and  then  compute  the  matrix  x  to  maximize  the  objective 
function  value  by  solving  the  following  maximization  problem: 

max  P(x),  (25) 

x(fc)ec 
1  <k<K 
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where  x(fc)  =  [xlk,x2k% *  •  •  ,xNk]T  and  C  =  {(xi,x2,  ■  ■  ■  ,xn),Xi  >  0,  i  =  1,  •  •  •  ,  N,  Y.iL\  xi  =  1} 
is  a  convex  set  in  KA\  This  solution  matrix  x  determines  the  search  results  according  to  the  numerical 
values  in  the  entries  of  x.  For  document  i,  the  document  j  is  most  similar  to  document  i  if  the  distance 
between  vectors  ,  k  =  1,  •  •  •  ,  K  and  Xjk,  k  =  1,  •  •  •  ,  K  is  the  smallest. 

Since  P(x)  is  a  quadratic  function  of  x  and  hence  continuous,  P(x)  achieves  its  maximum  within  the 
bounded  set  C.  It  is  easy  to  see 

Lemma  5.1  Let  M  =  maxi <ij<n  Wij  be  the  maximum  of  the  entries  of  the  global  similarity  matrix. 
Then  -P(x)  <  MK  for  all  x(fc)  G  C,  k  =  1,  •  ■  •  ,  K. 

Proof.  Since  Ya=  t  xi  =  1  and  xi  -2  0  f°r  all  i  =  1,  •  •  •  ,  N,  we  have 

K  K  I<  n  n 

pw  =  E  E  3'ik'Wij  ^  EE  MxiXj  =  M  EE  xjkxik  —  MK. 

k=  1  ij  k=  1  ij  k= 1  i=l  j= 1 


However,  P(x)  may  achieve  its  maximum  at  several  vectors.  For  example,  letting  K  =  1,  N  =  3  and 
W  be  the  identity  matrix  of  size  3  x  3,  we  know  the  maximum  is  1  when  x  =  (1, 0, 0),  (0, 1,  0),  (0, 0, 1). 
Certainly,  this  is  just  a  pathological  example  as  no  one  uses  an  identity  matrix  for  the  global  similarity 
matrix.  The  algorithm  Yahoo  uses  produces  a  unique  maximizer  each  time  as  it  will  be  discussed  below. 

To  solve  the  maximization  problem,  Yahoo  uses  the  Growth  Transformation  algorithm! — (GTA).  The 
most  important  advantages  of  this  algorithm  arc  that  it  guarantees  a  monotonic  increase  sequence  of 
the  objective  function  values  and  the  computation  is  stable  as  each  step  is  a  convex  combination  of  the 
previous  step.  For  simplicity,  we  assume  that  K  =  1.  In  fact,  Yahoo  engine  does  a  loop  to  go  through 
each  category  using  the  following  GTA. 


Algorithm  5.1  (the  Growth  Transformation)  Starting  with  a  trained  vector  x1 1 '  in  C,fork  =  1,2,--- 
one  computes 


x(fc+1) 

xj 


xfgf-FfxW) 


N 


EP>!t(xA 

i=  1 


(26) 


forj  =  !,-••  ,  N . 


This  algorithm  is  mainly  based  on  an  inequality  proved  in  [5].  For  convenience,  we  present  a  short 
version  of  the  proof  due  to  the  degree  d  =  2. 


Lemma  5.2  (Baum  and  Eagon,  1967)  Let  P(x)  be  a  homogeneous  polynomial  of  degree  d  with  nonneg¬ 
ative  coefficients.  Let  x  G  C  the  convex  defined  above.  For  any  x  G  C,  let 


y  =  , Vn)  with yj  = 


x3lLPW 


N 


d 


-,j  = 


.  N. 


i=  1 


Then  P(y)  >  P(x)  unless  y  =  x. 
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Proof.  For  convenience,  let  us  use  P(x)  =  xTICx  to  prove  this  inequality,  i.e.,  we  assume  that  d  =  2. 
Clearly,  P  is  a  homogeneous  polynomial  of  degree  d  =  2.  Since  the  entries  of  W  are  nonnegative,  P(x) 
satisfies  the  assumptions  of  this  lemma.  Note  that 


d_ 

dxi 


P(x)  =  2 


WijXj 


and 

N  o 

=  2P(x)‘ 

i=l  1 

Thus,  we  first  use  Cauchy-Schwarz’s  inequality  and  then  the  geometric  and  arithmetic  mean  inequality  to 
have 


P(X)  =  52  WijXiXj  =  J2  (WijViyj)1'3  (wTXiX^7~vJ3 

i,j=l,-,n  i,j= 1  '  \yiVj) 

'  n  \  /  n 

52  wayiVi )  (  t  ,rux>xi  ( zrzf] 

,,  I  /  \  i,  7=1  V'  !li!l  i  J 


< 


< 


1/3 


52  wayiyj 


\ij=  1 


El  ( Xi  Xj 
WijXiXj- - h  — 

,y=,  w 


2/3 


Since  m  =  xy  ^  WijXj/0(x),  we  have 


1=1 


1  f  Xi  Xj  \  1  P(x)  P(x) 

>  WijXiXj-  —  +  —  =  -  >  WijXiXjl^n - +  ™ - 

2  V2/i  vjJ  -  ;y  Ej=imjXj  E*=1  wvxi 


1= 1  *=1 


y 


s/nn  n  n  n  n  \ 

2  «W/(£  WikXk^  +  X]  X1  X]  WijXi/ (52  WkjXk)  ) 

V*=i  i=1 

\*=i  i=i  / 


=  P(x)- 


since  x£C.  Combining  the  above  discussion  together,  we  have 

p(x)  -  f  XI  wnviyj j  (p(x))2/3- 


n 

It  thus  follows  that  P(x)  <  Wijmyj  =  P(y).  Note  that  from  the  above  discussion,  the  inequalities 

*,i=t 

become  equalities  for  Cauchy-Schwarz’s  inequality  and  geometric  and  arithmetic  mean  inequality  when 
y  =  x.  Therefore,  we  have  an  strictly  inequality  when  y/x.  We  have  thus  completed  the  proof.  ■ 
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With  this  lemma,  we  arc  able  to  show  the  convergence  of  Algorithm  5.1.  Due  to  the  monotonically 
increasing  property  of  P(x.(k'1)  and  the  boundedness  of  P(x)  by  Lemma  5.1,  we  see  that  P(x®)  is 
convergent  to  P(x*)  for  some  point  x*.  In  practice,  the  iterations  may  end  in  a  finite  steps.  Note  that  for 
x  G  C,  VP(x)  f  0  and  each  component  of  VP(x)  is  strictly  positive.  Thus,  objective  function  values 
arc  increasing  until  P  achieves  its  maximum.  In  this  case,  by  Lemma  5.2,  we  have  x(n+ 1  =  xkn>  for 
some  integer  n. 

Also,  x^k> ,  k  >  1  arc  all  in  C  and  hence,  arc  bounded.  It  follows  that  there  exists  a  convergent 
subsequence.  Due  to  the  nonuniqueness  of  the  maximizers,  we  can  not  conclude  the  convergence  of  the 
whole  sequence. 

Suppose  that  W  is  a  positive  definite,  P(x)  is  a  convex  function,  when  x  and  y  defined  in  Lemma  5.2 
are  close  enough,  we  have 

VP(x)  •  (y  -  x)  >  0  (27) 

since  P(y)  >  P(x).  It  follows  that 

P(y)  -  P(x)  =  P(y  -  x)  +  2(y  -  x)Tlfx  =  P(y  -  x)  +  (y  -  x)TVP(x)  >  An||y  -  x||2,  (28) 

where  An  denotes  the  smallest  eigenvalue  of  W  and  ||z||  denotes  the  (2  norm  of  vector  z  e  Mn.  Therefore 
we  have 


Theorem  5.1  Suppose  that  W  is  symmetric  and  positive  definite.  Let  M  be  the  largest  value  of  the  entries 
ofW.  Let  An  be  the  smallest  eigenvalue  ofW.  Let  ySk\  k  >  I  be  the  sequence  from  Algorithm  5.1.  Then 
there  exists  a  convergent  subsequence,  say  x.(mk\  k  >  1  which  converge  to  a  maximizer  x*  and  satisfy 

OO 

J2  Hx^+i)  -  x^f  <  C  <  00, 
k= 1 


where  C  is  a  positive  constant  dependent  on  An  and  M. 

Proof.  We  have  discussed  that  the  sequence  x(k>  has  a  convergent  subsequence.  For  convenience,  let  us 
say  x(ki  converges  to  x*.  Thus,  x(/,"+ 1  -1  —  x:  k>  will  be  very  close  when  k  >  I\ 0  for  a  positive  inteeger  Kq. 
Thus, 

P(x(fc+L)  -  P(x(fc))  >  P(x(fc+1)  -  x^)  >  An||x(fc+1)  -  xW||2 

by  using  (28)  and  (27).  Hence, 


An  Hx(fc+1)  -  xWl|2  ^  M  ~  p(x(i"0))  <  M  <  OO. 

k>K0 


This  completes  the  proof.  ■ 

Although  Yahoo  does  not  discuss  much  mathematics  behind  their  generating  the  global  similarity 
matrix  W,  one  approach  called  matrix  completion  can  be  used.  The  matrix  completion  is  a  research  topic 
in  mathematics  which  is  recently  actively  studied.  It  stalls  from  the  well-known  Netflix  problem.  Indeed, 
Netflix  (cf.  [36])  made  available  publicly  such  a  set  of  data  with  about  105  known  entries  of  its  movie 
rating  matrix  of  size  about  5  x  105  times  2  x  105  and  challenged  the  research  community  to  complete 
the  movie  recommending  matrix  with  root  mean  square  error  less  than  0.8563.  This  matrix  completion 
problem  has  been  studied  actively  ever  since.  We  refer  the  following  papers,  [26],  [41],  and  [29]  for 
in-depth  study  and  the  references  therein. 
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6  Remarks 


We  have  the  following  remarks  in  order: 

Remark  6.1  We  will  not  discuss  other  link-based  ranking  systems  for  webpages  such  as  the  HITS  algo¬ 
rithm  invented  by  Jon  Kleinberg  (used  by  Ask.com),  the  IBM  CLEVER  project,  IDD  Information  Sendees 
designed  by  Robin  Li  who  found  Baidu  in  China,  the  TrustRank  algorithm  and  etc.. 

Remark  6.2  Next  we  discuss  how  to  in  parallel  update  the  RageRank  simultaneously  for  multiple  web¬ 
pages  and  keywords  and  names  of  documents  at  a  time.  We  need  the  following 

Lemma  6.1  ( Woodbury  matrix  identity )  Suppose  that  A  is  invertible  and  C  is  also  invertible.  IfA+U CV 
is  invertible,  then 


(A  +  UCV)-1  =  A-1  -  A-'UiC-1  +  V A~lXJ)~lV A-1 . 

The  above  formula  can  be  rewritten  in  the  following  form 

Lemma  6.2  Suppose  that  U  and  V  are  unitary  matrices  and  C  is  a  diagonal  matrix  containing  nonnega¬ 
tive  values  as  entries.  If  A +  U CV  is  invertible,  then 

(A  +  UCVT)-1  =  A~x  -  (VC^UT A  +  I)-1  A-1,  (29) 

where  CL  is  the  pseudo-inverse  ofC. 

Proof.  First  of  all,  we  have 

(. A  +  UCVT)~1  =  A'1  -  J4_1?7(C't  +  1 7r  A~1U)-1Vr  A-1 


by  using  the  same  direct  proof  for  Woodbury  formula.  Furthermore,  we  use  the  invertibility  ofU,  V  and 

A  to  rewrite  the  right-hand  side  of  the  above  formula  to  get  the  equality  in  (29).  ■ 

k 

For  the  Google  matrix  A  =  I  —  dM  and  UCVT  =  —d  XjvJ,  we  know  most  entries  of  matrix 

3= 1 
k 

XjvJ  of  size  N  x  N  are  zero.  There  is  a  nonzero  block  which  is  of  size  rtk  x  rik  with  integer 

3= 1 

rik  «  N  if  k  is  not  very  big.  For  convenience,  let  us  say  the  principal  block  of  size  rik  X  rik  are 

k 

nonzero.  It  is  easy  to  find  the  singular  value  decomposition  of  this  block,  -y>vj  =  U\C\V^ .  Then 

3= 1 

let  U  =  diag(U\,  J/v_nfc),  C  =  diag(C\ ,  0 N-nk)  and  V  =  diag(V i,  /jv-nfc).  where  lN-nk  is  the  identity 
matrix  of  size  N  —  n &  and  0 N-nk  is  zero  block  of  size  N  —  nj~.  In  general  C\  is  not  invertible.  We  can 
write  it  in  the  following  form:  if  C\  is  of  rank  rk  <  rik, 


Ci 


c  11  o' 

0  0  7 
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where  C\  \  is  diagonal  matrix  containing  all  the  nonzero  singular  values  of  —  XjvJ. 

3= 1 

From  (29),  we  need  to  find  (VC^UT  A  +  I)~l.  First,  we  have 

VC'UTA=\Vl  °  1  f^ii1  0  IK  0  }\lnk-dMu  —dM \2 


0  lN-nk  o  On—ti,  o  I N—nk.  —dM2l  lN-nk  ~  dM22 


for  a  matrix  C\  which  is  of  size  n X  n^.  Then  we  have 


vrUrr  A  =  [K  0  1  \lnk-dMu  —dM \2 

-  [  0  Otv— nj  [  -dM2 1  IN-nk  -  dM22\  ’ 

where  the  Google  matrix  M  is  written  in  terms  of  the  same  division  of  blocks  as  VC'U  ' .  Furthermore, 

VC]UT  A  +  I  =  -  dMn)  +  Ink  -CiM12 

[  0  IN-nk 

_'\Ci  0  1  (\dCfl  -dMu  -dM12\  \ 

0  lN-nk  V  0  0  N-nk  ) 


Finally,  we  have 


(ytfUT  A  +  I)-1  =  (l  +  d\Cl-  Mu  ~Ml2] )  1  K  r  0 
V  L  0  ®N-nk\)  [  0  In-i 


It  follows  the  updated  PageRank  vector  is 


P  +  =  (I-  dM  -  dUCVT)~ 

=  (I-dM)-1^l-  (i  +  d 


C{ -Mn  -Mi21\  1  \dC\  0 


0  On-t 


0  In-i 


(I  -  dM  f 


=  P-  [I  +  d 


C\  —  Mu  ~Ml2]\  1  \dC\  0 


0  0/v— r 


0  In-i: 


We  are  now  ready  to  present  an  algorithm  for  multiple  updates. 

Algorithm  6.1  (Multiple  Updates )  Choose  an  integer  K  large  enough  such  that  dh  is  within  a  given 
tolerance.  We  first  compute  U\,C\,V\  as  explained  above.  Then  compute  C\.  For  the  current  PageRank 
vector  P,  modify  it  to  be 

p  =  dC\  0  p 

0  lN-nk\ 


and  then  compute  the  following  iterations 


P+«P  -JV 

3=0 


j\c\-Mu  -Mi2Y 

0  0  N-nk_ 
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Similarly  we  can  discuss  how  to  compute  M+  efficiently  when  adding  a  new  webpage/document/link 
to  the  database.  We  omit  the  detail  here. 

Remark  6.3  The  importance  of  webpages  related  to  each  keyword,  document  or  webpage  can  also  be 
calculated  based  on  the  number  of  hits  from  the  world.  However,  this  is  very  easily  be  scrolled  up  by  some 
artifical  hits. 

Aknowledgment:  The  authors  would  like  to  thank  Dr.  John  Lavery  for  his  encouragement  and  support 
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